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Abstract 

In an early approach, we proposed a kinetic model with multiple translational temperature [K. Xu, 
H. Liu and J. Jiang, Phys. Fluids 19, 016101 (2007)], to simulate non-equilibrium flows. In this 
paper, instead of using three temperatures in x— , y—, and z-directions, we are going to further 
define the translational temperature as a second-order symmetric tensor. Based on a multiple stage 
BGK-type collision model and the Chapman-Enskog expansion, the corresponding macroscopic gas 
dynamics equations in three-dimensional space will be derived. The zeroth-order expansion gives 
the 10 moment closure equations of Levermore [CD. Levermore, J. Stat. Phys 83, pp.1021 (1996)]. 
To the lst-order expansion, the derived gas dynamic equations can be considered as a regularization 
of Levermore's 10 moments equations. The new gas dynamic equations have the same structure 
as the Navier-Stokes equations, but the stress strain relationship in the Navier-Stokes equations is 
replaced by an algebraic equation with temperature differences. At the same time, the heat flux, 
which is absent in Levermore's 10 moment closure, is recovered. As a result, both the viscous 
and the heat conduction terms are unified under a single anisotropic temperature concept. In the 
continuum flow regime, the new gas dynamic equations automatically recover the standard Navier- 
Stokes equations. The current gas dynamic equations are natural extension of the Navier-Stokes 
equations to the near continuum flow regime and can be used for microflow computations. Two 
examples, the force-driven Poiseuille flow and the Couette flow in the transition flow regime, are 
used to validate the model. Both analytical and numerical results are encouraging. 



1 Introduction 

The transport phenomena, i.e., mass, heat, and momentum transfer, in the different flow regime 
is of a great scientific and practical interest. The classification of various flow regimes is based 
on the dimensionless parameter, i.e., the Knudsen number, which is a measure of the degree of 
rarefaction of the medium. The Knudsen number Kn is defined as the ratio of the mean free path 
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to a characteristic length scale of the system. In the continuum flow regime where Kn < 0.001, the 
Navier-Stokes equations with linear relations between stress and strain and the Fourier's law for 
heat conduction are adequate to model the fluid behavior. For flows in the continuum-transition 
regime (0.1 < Kn < 1), the Navier-Stokes (NS) equations are known to be inadequate. This 
regime is important for many practical engineering problems. Hence, there is a strong desire and 
requirement for accurate models which give reliable solutions with low computational costs. 

One of the alternative approach to simulate the non-equilibrium flow is those based on the 
moment closures. Grad's 13 moment equations are one of the most important ones, which provide 
the time evolution of the non-equilibrium quantities, such as the stress and the heat flux [8]. 
However, due to its hyperbolic nature, these equations lead to a well known sub-shocks problem 
inside a shock layer as the Mach number is larger than a critical value. In order to improve 
the validity of the 13 moment equations, based on the Chapman-Enskog expansion, Struchtrup 
and Torrilhon introduced terms of super-Burnett order to the balance of pressure deviator and 
heat flux vector in the moment equations, and got the regularized 13 moment (R13) equations 
which have much better performance in the non-equilibrium flow regime [20]. Another well-known 
moment system is Levermore's 10 moment closure, which follows his hierarchy of non-perturbative 
moment closures with many desirable mathematical properties |13| . For example, these equations 
don't suffer from the closure-breakdown deficiencies, and they always give physically realizable 
solution due to non-negative gas distribution function. However, the 10 moment Gaussian closure 
has no heat flux even though it is proved that Navier-Stokes viscous terms can be recovered in 
the continuum flow regime. In an effort to extend the Gaussian closure to include higher-order 
effects, Groth et. al. formulated perturbative variants to the original moment closure with new 
extended fluid dynamic model |10j . And the most well studied of these closures is a 35- moment 
closure. Recently, McDonald and Groth took a Chapman-Enskog-type expansion about either the 
moment equations or the kinetic equation, and introduced the heat flux into Levermore's 10 moment 
closures and obtained extended fluid dynamic equations for non-equilibrium flow simulation [17] . 
The new system present improved results in the transition flow regime where the heat transfer has 
a significant effect. 

In recent years, we have concentrated on the development of numerical schemes for the near 
continuum flow simulation. In order to capture the non-equilibrium physics in the transitional flow 
regime, we have extended the gas-kinetic Navier-Stokes flow solver with the following developments 
[23]. First, a closed solution of the gas distribution function up to the NS order has been used to 
derive a generalized particle collision time, subsequently to obtain the extended viscosity and heat 
conduction coefficients [23]. Later, in order to describe the non-equilibrium flow related to the 
molecular rotational and vibrational degree of freedom, a multiple time relaxation kinetic scheme 
has been introduced for the shock structure calculations [28, 5j. Recently, the gas-kinetic scheme has 
been further extended to study the multiple translational temperature non-equilibrium [30]. The 
schemes developed in the above study give reasonable results in the transitional flow regime, such 
as the capturing of shock structure at different Mach numbers and flow phenomena which cannot be 
described properly by the NS equations. In the above simulations, the underlying physical model 
is a generalized BGK (GBGK) model [30], where the multiple stage relaxation processes have been 
considered. 

In the gas- kinetic schemes, the solutions are obtained without knowing the explicit macroscopic 
governing equations. In this paper, following the numerical procedure we are going to fill up the 
gap and derive the underlying macroscopic governing equations for the monatomic gas. In an 
early approach [30], we introduced the multiple translational temperature into the kinetic model, 
where the energy exchanges between x-, y-, and z-directions are modeled through the particle col- 
lision. Based on the above kinetic model, in one dimensional space the generalized NS equations 
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are derived, where the viscous term in the NS equations is replaced by the temperature relaxation 
term. In this paper, we will further develop such a model and regard the temperature as a sec- 
ond order tensor. Since the gas flow may settle to an equilibrium state through multiple stages 
[26], one of the reasonable assumption is to use a Gaussian distribution with multiple temper- 
ature as a middle state. Physically, this state corresponds to a gas with different temperature 
in different directions. Therefore, the thermal energy or particle random motion is represented 
through a symmetric temperature tensor 1^-. Due to historical reason, the temperature is de- 
fined as a thermodynamical variable, where there is no bulk fluid velocity variations. So, the 
temperature becomes a scalar concept. However, in the transition flow regime the gas has large 
bulk velocity variation, and there is no enough particle collisions to equalize the random particle 
motion. In order to construct gas dynamic equations, it is justified to extend the temperature 
from a scalar concept to a tensor. The commutable property of the random particle velocity de- 
termines the temperature to be a symmetric tensor. Actually, this kind of non-equilibrium gas 
property has been routinely extracted from the DSMC solutions. So, based on the physical model 
with the Gaussian distribution as a middle state between the real gas distribution function / and 
the equilibrium Maxwellian f eq and the strategy used in the construction of the kinetic scheme, 
we are going to derive the corresponding macroscopic governing equations. Surprisingly, the ob- 
tained gas dynamic equations become regularization of Levermore's 10 moment Gaussian closure, 
where additional viscous and heat conduction terms are obtained. The structure of gas dynamic 
equations are almost identical to the Navier-Stokes equations, but the constitutive relationship 
of the NS one, i.e., aij = —pRT eq 5ij + p(diUj + djUi — ^dkUkSij), is replaced by the new one 
<Jij = —pRTij + pR(T^ — Tij). At the same time, the heat flux depends on the gradient of the 
temperature T^. 

This paper is organized in the following. Section 2 is about the introduction of kinetic equation 
and the generalized particle collision model. At the same time, the 10 moment closure and the 
generalized gas dynamic equations will be presented. Section 3 is about the applications of the new 
gas dynamic equations to two flow problems in the near continuum regime. The last section is the 
conclusion. 

2 Generalized Gas Dynamic Equations 

In this section, we first review the particle collision model, introduce the Gaussian closure, and 
derive the new gas dynamic equations. A monatomic gas will be considered in this paper. 

2.1 Two stage gas-kinetic collision model 

The gas-kinetic Bhatnagar-Gross-Krook (BGK) model has the form [3], 



where the particle distribution function / is a function of time t, spatial location Xi, and particle 
velocity tij. The left hand side of the above equation represents the free streaming of molecules in 
space, and the right side denotes the simplified collision term of the Boltzmann equation. In the 
BGK model, the collision operator involves a single relaxation time r for a non-equilibrium state 
to evolve to an equilibrium one f eq , which is an isotropic Gaussian 



(1) 



d t f + u l d l f=(f eq -f)/r, 
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where p is the density, T eq is the equilibrium temperature, and U% is the averaged macroscopic 
fluid velocity. Traditionally, based on the above BGK model, the Navier-Stokes and higher-order 
equations, such as Burnett and Super-Burnett, can be derived [3CE]. Unfortunately, these higher- 
order equations have intrinsic physical and mathematical problems in the transitional flow regime. 
In general, the BGK collision term is valid only for flows close to the thermal equilibrium one. In 
order to extend the capacity of the BGK model to the non-equilibrium flow regime, we can re-write 
the collision term of the BGK model into two physical sub-processes, 

(2) d t f + uAf = (g- f)/r + (f eq - g)/r, 

where g is a middle state between / and f eq , see figure 1. The DSMC solutions show that the 
randomness of particle motion depends on the spatial direction. So, a nature assumption about the 
middle state is a state with multiple temperature. In the above equation, the term (f eq — g)/r has 
no direct connection with /, therefore, we can consider it as a source term in the above generalized 
BGK (GBGK) model, 

(3) d t f + utdif = (g- f)/r + Q, 

where Q = (f eq — g) jr for the monatomic gas. 

In an early approach, we have assumed that the middle state g has individual temperature in 
x— , y— and z— directions, and it has the following form in a two dimensional case, 

n v 1/2 f] x 1/2 /, x 1/2 

(4) 9 = P -) F) (-] exp \-l x { u -Uf-l y {v-Vf-l z w< 



7T / \ 7T / \ 7T 

Here l x = 1/(2RT X ), l y = l/(2RT y ), and l z = \/{2RT z ) are related to the translational temperature 
T x ,T y , and T z in x, y, and z directions. Here R is the gas constant. 

In order to solve the above kinetic equation numerically to determine the time evolution of 
the macroscopic physical quantities (p,U,V,T x ,T y ,T z ), we proposed the following scheme. First, 
we expand the gas distribution function around the multiple temperature state g, such as / = 
g — r{gt + Uidig), from which the numerical fluxes across each cell interface are evaluated. Second, 
the source term Q is integrated in a time step inside each computational cell and is regarded as 
a source term. The physical effect of the source term is to equalize the temperature in different 
directions. The above processes are composed of the relaxations from the non-equilibrium state / 
to a multiple temperature state g, then g converges to an equal temperature Maxwellian. These 
two processes may have different relaxation time scales, and the specific formulation of g depends 
on the flow problems. According to the above model, we have derived the hydrodynamic equations 
in one dimensional space [30], where the stress-strain relationship in the Navier-Stokes equations 
is replaced by the temperature relaxation. In the continuum flow regime, where the middle state 
is close to the Maxwellian, the standard Navier-Stokes equations can be recovered. The numerical 
tests presented in |30| verified the validity of the above model. In this paper, we are going to 
use multiple temperature Gaussian as the middle state, and according to the similar procedure to 
derive general gas dynamic equations in three dimensional space. 

To approximate the Boltzmann particle collision term by multiple BGK-type sub-processes 
have been investigated before by many authors, such as Callaway [6j, Gorban and Karlin [9], and 
Levermore |13j . For the DSMC method, to include the rotation and vibration modes has been done 
through the BGK-type relaxation model as well. As realized in [13], for Maxwellian molecules, 
the use of the multiple BGK-type collision term is correct even if the flow is far away from the 
equilibrium. This suggests that the generalized BGK operator may be a legitimate approximation 
to the collision term for use in the near continuum flow regime. 
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2.2 Generalized gas dynamic equations: zeroth order 

For the rarefied flow simulation, a generalized middle state g between / and f eq can be a Gaussian 
distribution. The Gaussian distribution appears to have been derived in the early work by Maxwell 
[16j . and then re-discovered by many researchers, such as Holway [12J and Levermore [13J. In 
this paper, we directly define the gas temperature as a tensor instead of a scalar. The traditional 
temperature concept is coming from thermodynamics, where the local bulk fluid velocity deviation 
is absent. However, for the transport equations, the temperature basically represents the degree 
of randomness of particle motion. In order to dynamically equalize the temperature, there should 
have enough particle collision. In the near continuum flow regime, the number of particle collision 
inside a microscale device, is limited. It is most likely that the particles will keep the non-isotropic 
particle randomness. 

One of the middle state we can use between / and f eq is the Gaussian distribution, 

9 = , 9 = exp(-^ - U^RTij)- 1 ^ - Uj)), 
/det(27ri?T ii ) 2 



where Ty is the positive definite temperature matrix, which is related to the thermal energy of the 
particle motion pRTij = J(ui — Ui){uj — Uj)gdu. Due to the commutable property between particle 
randomness velocities, such as (itj — Ui)(uj — Uj) = (uj — Uj)(ui — U), must be a symmetric 
tensor. In the Levermore's approach [T3], he developed a non-perturbative method, where the gas 
distribution function / is assumed to be equal to g. Based on the kinetic equation with the 
assumption f = g, i.e., 

d t g + Uidig = Q, 

taking the moments 

tp = (l,Ui,UiUj) T , 

on the above equation gives the following Gaussian closure, 

(5) d tP + d k (pU k ) = o, 

(6) d t (pU) + d k [p{UU k + RT lk )} = 0, 

(7) dMUUj + RTij)} + d k [p(UUjU k + RU k Tij + RU t T jk + RU 3 T kl )} = ^ P R(T eq 5 lJ - Ty). 

These equations are actually the zeroth-order Chapman-Enskog expansion, i.e., / = g, for the 
generalized BGK model. In the above equations, the source term on the right hand side is coming 
from the term Q in the generalized BGK model. The equilibrium temperature T eq is defined by 

T eq = ^Tr(T^). 

If the state g is equal to the Maxwellian, i.e., g = f eq , the above equations reduce to the Euler 
equations. Based on the above equations, Levermore and Morokoff shows that if the initial data of 
T{j is symmetric positive definite, then it remains so [14J. 

Even without heat conduction terms, Levermore showed that the above equations recover the 
Navier-Stokes viscous terms in the continuum flow regime |13| . For the above system ([5])-([7|). 
the left hand side equations have a complete real eigenvalues and eigenvectors [2J [TT], and the 
system is strictly hyperbolic. The temperature difference in different direction, such as the T{j, 
basically shows that the sound speed depends on the spatial direction, the so-called anisotropic 
wave propagation due to the non-isotropic gas property. 
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In order to further introduce heat flux into the above 10 moment closure, McDonald and Groth 
recently took a Chapman-Enskog-type expansion about the Gaussian moment equations [17J. The 
heat conduction is added for the thermal temperature equation, 

(8) d t T i:j + 8 k {U k T %J ) + T lk d k Uj + T jk d k Ui + 8 k Q ijk = ~(T e %- - Tij), 
where the heat flux Q has the form 

T 

Qijk = —^[TkidiTij + TjidiT ik + TftdiTjk]. 
2.3 Generalized gas dynamic equations: first order 

The generalized BGK model includes two relaxation process. One is from / to the Gaussian g, 
and the other is from g to an equilibrium state f eq . In the last section, the distribution function 
/ is assumed to be equal to g, and the 10 moment closure equations are derived. However, as 
presented in figure 1, the real distribution function should be different from g, and the process 
of relaxation from f to g has to be considered. In the past years, we have developed gas-kinetic 
schemes based on the generalized BGK model, where a gas distribution function / around g has 
been constructed and used to evaluate the numerical fluxes in a finite volume scheme |30| . The 
schemes present reasonable numerical solutions in the near continuum flow regime, such as the 
micro-channel flow computation [26J. In the following, we are going to derive the corresponding 
macroscopic governing equations underlying the gas-kinetic scheme. The method used here can be 
regarded as the Chapman-Enskog expansion or iterative expansion [19], which are equivalent to 
each other. 

The solution / around the Gaussian g is constructed using the iterative expansion to the 1st- 
order order [19] . 

(9) f = g~ r{d t g + uAg) + tQ, 

where in the kinetic scheme dtg is determined using the compatibility condition 

(10) J ip(d t g + Uidig)du = J ipQdu, 

which is exactly the 10 moment closure ([5])- ([7]). Substituting the distribution function / in ([9]) into 
the BGK model ([3]), the equation becomes 

(11) dtg + Uidig = r(d 2 t g + 2uid t dig + UiUjdidjg) + Q- r(Q t Q + UidiQ). 



Taking the moments ip to the above equation and using Eq. (ll0p to express the time derivative in 
terms of the spatial derivative, we can get the following macroscopic equations, 

(12) d t p + d k ( P u k ) = o, 

(13) d t (pUi) + d k [p{UiU k + RT ik )] = d k [ P R(T e(1 5 kl - T ki )], 
dMUiUj + RTij)] + dklpiUiUjUk + RU k Tij + RUiT, jk + RUjT H )] 

+d k {pR[U k (T eq 8 l3 - T^) + U t (T e % k - T jk ) + U 3 {T e H ki - T kl )}} 

(14) +d k { T pR 2 ^-{T kl d l T l , J + T a diT jk + TfldiTu)}. 
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The above equations have been written in a similar way as the Navier-Stokes equations. The differ- 
ences between the above equations and the 10 moment closure are the additional terms appeared 
on the right hand sides of the momentum and energy equations. It is interesting to see that the 
corresponding heat conduction term derived above has the same form as that obtained by Mc- 
Donald and Groth, even though they are obtained through different considerations [T7j. Based on 
the BGK-type collision model, a unit Prandtl number is obtained for the heat conduction term. 
However, since we believe that up to the NS order the structure of the gas dynamic equations 
will not be changed due to the BGK collision term or the exact Boltzmann collision model, we 
add the Prandtl number in the above corresponding heat flux term. Since the viscosity and heat 
conduction coefficients are the concepts for the continuum flow, the particle collision time r in the 
above equations is defined according to the result in the continuum regime, 

r = p/(pRT eq ), 

where p is the dynamical viscosity coefficient. Certainly, in the rarefied flow regime, the correspond- 
ing viscosity coefficient has to be modified [23]. Also, attention should be paid on the relaxation 
term in the energy equation, where additional 2 appears. 

The above equations can be re-arranged to get the time evolution equation for the thermal 
energy pRTij, 

dtWij] + d k \pU k Tij] 
= -p(T^d ki -T ij ) 

T 

-p[Tkjd k Ui + T ki d k Uj] (a) 
+d k [pU k {T eq 5 ij - T i:j )} + p{T eq 5 jk - T jk )d k Ui + p{T eq 5 k% - T ki )d k U 3 (b) 

(15) +d k {rp^(T kl d l T ij + T il d l T jk + T jl d l T ki )} (c). 

In comparison with the Navier-Stokes equations, the right hand side of the equation (|15p has a 
clear physical meaning. Here (a), the forces multiplied by fluid deformation, which is the heating 
and cooling of the fluid by compression or expansion and this term represents a reversible process, 
(b). viscous dissipation, which is responsible for heat generation, and it is always positive and 
produces internal energy. This is an irreversible process, (c). heat conduction term, which is also 
irreversible. 

The relaxation parameter r controls the distance between the non-equilibrium state / and the 
equilibrium one f eq . The current method for the derivation of the gas dynamic equations can be 
also used to other system, such as these with a more complicated non-equilibrium middle state g, 
such as the distributions with 14 or 26 moments. However, a distribution function g with higher- 
order terms may correspond to macroscopic governing equations without clear physical meanings 
for the higher-order terms. 

The generalized gas dynamic equations (|12p - (114p have the same left hand side as the 10 moment 
closure equations ©-(0) [14J. And the fluxes on the left have a complete eigenvalues and eigen- 
vectors [4j [TT| I18j . In other words, the left hand side is the hyperbolic part. For the right hand 
side, besides the heat flux as recently derived by McDonald and Groth [17], additional dissipative 
terms, i.e., pR(T eg 5ij — T^) in the momentum equations, and pR[U k (T eq 5ij — T^) + Ui(T eg 5j k — 
Tj k ) + Uj(T eq 5 k i — T k i)] in the thermal energy equation, have been obtained. Furthermore, in the 
present model the relaxation time appearing in the first term on the right hand side of Eq. (|14p is 
r/2, while the corresponding term in the model by McDonald and Groth is r. In our two-stage 
collision model. If we compare the above equations (|12p - (|14p with the Navier-Stokes equations, 
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where the energy equation can be also separated in a direction-by-direction componentwise form, 
we can immediately realize that the corresponding viscous term is given by 

Therefore, the constitutive relationship for the new gas dynamic equations is 

(16) ay = —pRTij + pR(T eq 5 i3 - Tij ) . 

In the following, we are going to show that the new constitutive relationship can recover the 
standard NS formulation in the continuum flow regime. 

In the continuum flow regime, the Gaussian distribution g will come to the same state as the 
equilibrium one f eq , see figure 1. In this case, the lst-order expansion of / presented in this 
section will be expanded basically around f eq . To the leading order of small r, Eq.(|14j) gives the 
temperature deviation, 

pRiT^-Tij) « 

(17) (r/2){d t [piUiUj + RTij)} + d k \p{UiUjU k + RU k T %J + RU % T jk + RUjT ki }]}. 

In the above equation, when applying the equilibrium condition Tu = T eq and = T eq 5ij, and 
using the Euler equations to replace the temporal derivative by spatial derivative, we can get 

pR(T eq 5ij - T^) 

(18) « (T^pRT^&Uj + djUi - ^d k U k 5 l3 ] = {ii/2)[diVj + OjUi - ^d k U k 5 i3 ]. 

Therefore, the constitutive relationship (|16p becomes 

ay = - P RT eq 5 i] +2pR{T eq 5 i] -T lj ) 

2 

(19) = -pSij + p[diUj + djUi - -d k U k 5ij}, 

which is exactly the Navier-Stokes stress-strain relationship. In the Navier-Stokes equations, the 
stress becomes symmetric tensor which is constructed through rational mechanical analysis. How- 
ever, in the new gas dynamic equations, the "stress" tensor is automatically a symmetric tensor due 
to the commutable property of random particle velocities, i.e., (ui — Ui)(uj — Uj) = [uj — Uj)(ui—Uj). 
In order words, the present formulation gives a microscopic interpretation of the origin of symmetric 
stress tensor. At the same time, the heat flux transported in A;— direction for the thermal energy 
pRT^ becomes 

Qtej = p r (TkidiTij + TudiTjk + TjidiTki). 

These results are consistent with the early analysis for the one-dimensional flow [30J. 

In summary, the generalized gas dynamic equations have been derived based on the multiple 
stage BGK model. With the added Prandtl number Pr and the introduction of dynamic viscosity 
in the determination of relaxation parameter r, the gas dynamic equations derived in this section 
are a closed system, which is a natural extension of the Navier-Stokes equations. Theoretically, the 
new gas dynamic equations cover a wider flow regime than that of the Navier-Stokes equations. 
The new constitutive relationship now has a microscopic physical basis. The viscous term involves 
only first-order derivatives of the flow variables. The ability to treat non-equilibrium flow problems 
without evaluating higher than first order derivatives would prove very advantages numerically. The 
method should be relatively insensitive to irregularity of the grid, the straight-forward computing 
using Discontinuous Galerkin framework, and the easy implementation of boundary conditions. 
Due to the low order equations, the less communication between the computational cells (small 
stencil) makes its numerical scheme more efficient to implement on parallel computing architecture. 



8 



3 Solutions of the generalized gas dynamic equations 



In this section, we are going to present two flow problems in the near continuum flow regime for the 
validation of the generalized gas dynamic equations derived in the last section. Both test problems 
are microchannel flows, but with different flow speed and non-equilibrium properties. 



3.1 Analytic solution of force-driven Poiseuille flow in near continuum flow 
regime 

In this subsection, we will apply the generalized gas dynamic equations to rarefied gas flows, such 
as the case of force driven Poiseuille flow between two parallel plates. Both Direct Simulation 
Monte Carlo (DSMC) and kinetic theory have shown that even with a small Kn, the pressure and 
temperature profiles in this flow exhibit a different qualitative behavior from those predicted by 
the Navier-Stokes equations [151 Ell Ell El EU EH]. Therefore, this flow can serve as a good test 
problem for any extended hydrodynamic equations intended for non-continuum flow computation. 

In our test, the two walls of the channel locate at y = ±H/2, respectively, and the force 
a = (a, 0, 0) is along the x direction. The flow is assumed to be unidirectional, i.e, d x (p = for 
any variable <j), and the velocity has only an x-component in the laminar and stationary case, i.e., 
U m = (U, 0, 0). Under such conditions, the generalized gas dynamic equations (Tl2l) - (ri4ll reduce to 

2 ^k_ pa = , 
ay 



d_ 

dy 



— (2P W -P) = 0, 



2_±fo P P dd xx \ _ 2{P XX -P) dU 

Pr dy { xy dy + yy dy ) ~ r + ^ xy dy ' 
3 d d8 yy \ _ 2(P yy -P) 



and 



Pr dy \ TPyy dy 

p d6 zz \ _ 2{P ZZ -P) 
Pvdy{ yy dy)~ r 

1 d ( d0 xv d0 vv \ 2P XV , s dXJ 



where pj = p9ij is the pressure tensor, and P = p8 is the pressure, with Oij = RTij and 9 = RT eq . 

It is difficult to obtain an analytical solution from the above nonlinear system. Here we will try 
to find an approximate solution using a perturbation method similar to that used in Ref. [22]. To 
this end, we first introduce the following dimensionless variables: 



v n-A p. -Eii o - e Ji fr-JL +-L 



where the variables with subscript represent the corresponding reference quantities. With these 
dimensionless variables, the system can be rewritten as 

(20) *%BL - ep = 0, 

dy 

(21) 2P yy -P = C, 
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— 2fP xy —^ + fP yy —^ = 2Pi5 2 ^ + AYv8P xy — , 

dy\ y dy dy r dy 



(23) ± ( fP ™lf) = 2P ' S 



T 




•2 (Pyy ~~ 


P) 


T 




■2 (Pzz ~ 


P) 


T 



( 24 ) 40 ( ^W^f ) = 2P ^ : 
and 

(25) c# +TFxy dy ) " ^ r +C ™ dy' 
where C is a constant, and e and 5 are two dimensionless parameters given by 

a . H fW 1 
d 



2^0 ' rov 7 ^ V2Kn' 

Now we assume that the force acceleration a is small such that e <C 1, then we can expand the 
dimensionless flow quantities in powers of e (we will omit the hat for simplicity hereafter): 

/0 = l + e 2 p (2) +o(e 4 ); U = eU (i) + (e 3 ), C = l + e 2 CV\ 
P xy = eP$+0(e 3 ), e xy = eeS + 0(e 3 ), 
P u = l + e 2 P® + 0(e 4 ), 6 ti = l + e 2 e ( S ) + 0(e 4 ), for i = x,y,z, 



and 



P= l +e 2p(2) +0(e 4 )) 0=1+e 2^(2) +o(e 4 )) r = 1 + e 2 r (2) +0(e 4 



The odd or even properties of the variables as functions of e are based on their symmetric properties 
in terms of the acceleration a. In general, the odd(even) velocity moments of the distribution 
function / are also odd (even) functions of e. Furthermore, from the definition Py = pOij we can 
obtain some useful relation which will be used later: 

(26) p W=pW-flW = pW-(9« PW=0&\ for i = X ,y,z. 

Substituting the above expansions into the nondimensional system, we can obtain the first order 
differential equation in e. Eqs. (f20|) and (f25|) go to 



dP 



(i) 



1 = 0, 



Pr5^- = ^-2Pr«5 2 pW, 



dy 

dy dy 
which give 

(27) ' Pg) =0^ =y, 

(28) = -by 2 + C, 

where C is a constant. Here, we have made use of the symmetry property of the velocity profile 
about y = 0. 
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From Eqs. (|2Tj) ([22]) . ([23]) . and (|24p. we can obtain the second order equations in e: 
(29) 2PW_ P (2) =C (2) > 



d¥ 2 ) 



(30) = 2Pr<5 2 (0i - #( 2 )) - 1 - SPr^V , 



(31) 3 ^L = 2P^(0g)-^) 



.<^0 
dy 
and 

(32) = 2Pr<5 2 (flg) - 9^ 



where the results up to the first order Eqs. ([27]) and (|28j) have been used. Summing up the above 
three equations we have 

d 2 



(33) f- 2 {?>^ + 29$) = -1 - 8Pr5V . 



From Eqs. (j3"Tj) and (f33j) . we can get 



„ f2l 2 . ,{VW t .\ 2K 2 4 43 2 231 

6> (2j = — A cosh K v v 4 H v 2 H ^ + 

3 ^3 y y 15 y 5(T 125K 2 

/.W „ , {VW„ \ 2K 2 , 77 2 693 

0}f/ = ,4 cosh Ky y A y 2 » + B. 

yy \ 3 J 15 50 250K 2 

Since 9^ is obtained, the other two temperatures 9^x and 9iz from Eqs. ([30]) and ([32]) can be 
constructed, respectively. 

Now let's find the second order pressure 

p(2). From Eqs. ([261) and (|29|) . we can obtain 
p( 2 ) = + 2(0( 2 ) - eg) = - ™* cosh ( y-^Ky) + ^y 2 + 231 



yy ' 3^3 y J 5 " 25K 2 ' 

Then, p( 2 ) can be determined from (|26p as 

f2) n (2) />f2^ ^2) , fy/lO Tr \ 2K 2 4 197 2 924 

p m = P m _ ,( 2)) = c m _ _ cosh \_ Ky j + + _ y2 + __ _ B 

With the above results, we finally get the approximate solutions of the problem: 

(35) L = 1 + ^ P m > 

(36) ^ = l +e 2^(2). 

Jo 

With these results, we are able to make some discussions on the velocity, pressure, and temper- 
ature profiles for the force-driven Poiseuille flow problem. First, it is clear that the velocity profile 
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is parabolic. But for the pressure and temperature, their profiles may be complicated due to the 
presence of the hyperbolic cosine function. 

It is interesting to compare the approximate solutions of the present gas dynamic equations to 
those of the Navier-Stokes-Fourier (NSF) equations. For the same flow problem, the NSF equations 
reduce to 

d ( dU\ 



(40) 7m = - eS {W +u " p=constl r„ = 1+€ 



< 37 » Ty{'^) + " a = ' 

< 38 > f=°> 

d /, dU\ fdU\ 2 

Using the similar perturbation method, we can get the following approximate solutions for the NSF 
equations: 

U - ,xfy\ 2 lT T ^ = 1 + e 2^(2) 5 

where 0® f = -(2/15)Pr<5 2 j/ 4 + D, and D is a constant. It is shown that the velocity profile of the 
NSF equations is also parabolic, which is the same as the that of the generalized hydrodynamic 
equations. However, the NSF equations give a constant pressure, which is qualitatively different 

(2) 

from the predictions of the new proposed model. For the temperature, it is evident that 9 n ^ is 
only part of 9^ 2 \ and has one local maximum at y = 0. On the other hand, the presence of terms 
of cosh(y) and y 2 in 9^ of the new model makes it possible to have a local minimum at y = 0, due 
to [9 {2) }' y =o = and [0 (2) ]£=o = -K 2 (W0A + 216)/135 + 1.72 > as A < 2.322K~ 2 - 2.16. 

To see this more clearly, we present the pressure and temperature profiles for the case of Kn=0.1 
which was studied extensively using DSMC method [31J. The constants A and B in 9 and P of 
the present hydrodynamic model are obtained by enforcing the values of 9 at y = and —H/2 to 
be the same as the DSMC data, and then C is obtained by enforcing P(y = 0) to be the DSMC 
value. The NSF solutions are obtained by enforcing their pressure and temperature values to be 
identical to the DSMC data at y = 0. In Figj2j the reduced temperature and pressure variations, 
q(2) _ ( T / To _ x y e 2 and p (2) _ (p/p Q _ i)/ e 2 > are s hown for both the NSF equations and the 

present model. It is seen that the temperature and pressure profiles of the present model are in 
qualitatively agreement with the DSMC data. For example, the temperature takes a bimodal shape 
and exhibits a local minimum at y = 0, and the pressure has two local maximums near the two 
walls. These critical flow behaviors are absent in the profiles of the NSF equations. For example, 
the temperature minimum only appears on the super-Burnett order if the traditional BGK collision 
model is used to construct the gas dynamic equations 125] . These observations demonstrate the 
fundamental difference between the present hydrodynamic model and the NSF equations. 



3.2 Numerical solution for Couette Flows 

For the new gas dynamic equations, a corresponding gas-kinetic scheme can be developed [30]. In 
order to further test the validity of the new governing equations, we will study the planar Couette 
flow here pp. The height of the Couette system is ho = 50nm. The wall's temperatures are fixed 
at To = 273-fT and at equilibrium, the pressure of the gas is 1 atm. Given that the walls distance 
is less than a mean-free path and the relative wall speed is high, the gas system will be strongly 
out of equilibrium. Specifically, the velocity distribution for the particles is non-Maxwellian. Same 
as the case in [T], the temperature along the channel is defined as T x and the one perpendicular 
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to walls is T z . With the fixed Knudsen number Kn= 1.25, and various wall velocities, such as the 
Mach number M = 0.5, 1.0, and 1.5, the simulation results from the new gas dynamic equations are 
shown in FigEl where the DSMC solutions are also included [1J. Clearly, at high Knudsen number 
and Mach numbers, the temperature is anisotropic. The new gas dynamic equations basically are 
capable in capturing this kind of non-equilibrium flow phenomena. Also, it is fully necessary for 
any gas dynamic equations to consider the temperature as a tensor for the non-equilibrium system. 

4 Conclusion 

Based on the multiple stage particle collision BGK model and the Gaussian distribution function as 
the middle state, the generalized gas dynamic equations have been derived. Since the gas temper- 
ature basically represents the molecular random motion, the direct extension of the temperature 
concept from a scalar to a second-order symmetric tensor T^ is physically reasonable. In the 
non-equilibrium flow regime, the randomness of the particle distribution indeed depends on the 
spatial orientation. The new gas dynamic equations have the same structure as the Navier-Stokes 
equatons, but the NS constitutive relationship, 

a i& = -pBT^Sij + iM&Uj + djUi - -dkUkSij) 

is replaced by 

At the same time, the heat flux in the k— direction for the transport of thermal energy pRTij 
becomes 

1 ki i = p r (TkidiTij + TudiTjk + TjidiTki). 

In the continuum flow regime, the generalized constitutive relationship and the heat flux term go 
back to the corresponding Navier-Stokes formulations. The new gas dynamic equations can be 
regarded as a regularization of Levermore's 10 moment closure [13]. The gas dynamic equations 
have a wider applicable flow regime than that of the NS equations. They capture the time evolution 
of the anisotropic non-equilibrium flow variables, as demonstrated in our examples. 

The traditional temperature concept is coming from thermodynamics, where there is no anisotropic 
particle random motion in space. However, for the non-equilibrium flow transport, due to the inad- 
equate of particle collision the random molecule motion can become easily anisotropic. To directly 
consider the temperature as a tensor rather than a scalar is a reasonable description for the non- 
equilibrium flow. For a dilute gas, due to the lack of long range particle interaction, the randomness 
particle motion is the only source for the dissipation in the system. Under the new definition of the 
temperature, all dissipative effects in a dilute gas system, such as the viscosity and heat conduction, 
can be unified under the same concept Tij. The current gas dynamic equations can be useful in the 
study of microflows. The further development of the gas dynamic equations based on the kinetic 
model and its scheme in [27], which are valid for the compressible shock waves, will be conducted 
in the near future. 
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Appendix: Gas Dynamic Equations in 2D Space 



As a special application, we consider 2D case, where there are 7 unknowns (p, U, V, T xx ,T xy , T yy ,T zz ) 
The generalized gas dynamic equations for those unknowns are 

(41) 



d\J | dF | dG 

dt dx dy 



dx 



dy 



where U is the vector macroscopic flow variables, and F and G are x and y direction flux vectors 
given by 

/ P \ ( P u \ 

p(U 2 + RT XX ) 



U 



P 

P u 

pV 

p(U 2 + RT XX ) 
p(UV + RT xy ) 
p(V 2 + RT yy ) 
pRT zz 



p{UV + RT xy ) 
p(U 3 + WRT XX ) 
p{U 2 V + 2URT xy + VBT S 
p{UV 2 + URT yy + 2VRT, 
P URT ZZ 



xy) 



G 



E v = 



P V 

p(UV + RT xy ) 

P (V 2 + RTyy) 

p{U 2 V + VRT XX + 2URT xy ) 
p(UV 2 + URT y y + 2VRT xy ) 
p(V 3 + 3VRT yy ) 
P VRT ZZ 





pR(T e i5 xx 
pR{T e «5 xy 



l xy) 



P R[V(T e "5 xx - T x 



3pR[U(T q 6 xx T xx ) + tR(T xx 8 x T xx + T X ydyT xx )\ 

(T ^5 X y T X y -f" T R(T X yd X T XX ~\- TyyOyT XX ~\~ 2*T XX X T X y ~\~ 2<T X y Oy T X y ) ) ] 



pR\jJ (T ^3yy Tyy) H - 2V (T ^S X y T X y} ~\- T R(2<T X yO X T X y ~\- 2iTyydyT X y ~\~ T XX 3 X Tyy ~\~ T X yOyTyy^ 

zz T zz ) + t R(T xx d x T zz + T X ydyT zz )\ J 





P R{T e «6 xy 

pR(T e «5yy 



-xy) 



pR[V(T ^5 XX T xx ) + 2U (T e< ^8 xy T xy ) + t R(2T xy d y T xy + 2T xx d x T X y + T yy dyT xx + T X yd x T xx )\ 

pR\jJ (T ^3yy Tyy) ~t~ 2V (T ^S X y T X y ~\~ T R(T X y3yTyy ~\~ T XX X Tyy ~\~ 2Tyy3yT X y ~\~ 2T X yQ X T X y)^j\ 

3pR[V{T^5yy-Tyy)+TR{TyydyTyy+T X yd X Tyy)] 



V 



and the source term is 



&ZZ T ZZ ) + TR(TyyOyT ZZ + T X y X T ZZ ) ] 







2pR(T e i5 xx 
2pR(T e «5 xy 

2 P R{T^8yy 

\ 2pR(T e *6 ee 



) 



)/r 
)/r 

Uy)/T 



-xy 

L yy 



T zz )/t ) 
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In the above equations, the equilibrium temperature 

T e i = T £ i8 xx = T«>8 m = T^S XZ = ^tr^), 

and 
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V/H v/H 



Figure 2: Reduced temperature (up) and pressure (down) variations at Kn = 0.1, where the results 
from the DSMC [31], the current equations, and the Navier-Stokes equations are presented. 
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Figure 3: Ratio of the temperature components T x /T z vs the vertical position Z = z/ho in planar 
Couette flow for Knudsen number Kn=1.25, and Mach numbers M = 0.5, 1.0, and 1.5. The solid 
lines are the solutions of the new gas dynamic equations and the circles are the DSMC solutions 
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